Exact spectrum of the Lipkin-Meshkov-Glick model in the thermodynamic limit and 

finite-size corrections 
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The spectrum of the Lipkin-Meshkov-Glick model is exactly derived in the thermodynamic limit 
by means of a spin-coherent-state formalism. In the first step, a classical analysis allows one to 
distinguish between four distinct regions in the parameter space according to the nature of the 
singularities arising in the classical energy surface; these correspond to spectral critical points. The 
eigenfunctions are then analyzed more precisely in terms of the associated roots of the Majorana 
polynomial, leading to exact expressions for the density of states in the thermodynamic limit. Finite- 
size effects are also analyzed, leading in particular to logarithmic corrections near the singularities 
occurring in the spectrum. Finally, we also compute expectation values of the spin operators in 
a semiclassical analysis in order to illustrate some subtle effects occurring in one region of the 
parameter space. 
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I. INTRODUCTION 

The Lipkin-Meshkov-Glick (LMG) model was pro- 
posed in 1965 to describe shape phase transitions in 
nuclei 1 ' 2,5 . This model is often used to describe the 
magnetic properties of molecules such as Mni2 acetate 1 . 
However, it also captures the physics of interacting 
bosons in a double-well-like structure ' 6 and is thus rel- 
evant to (two-mode) Bose-Einstein condensates' as well 
as Josephson junctions. It has also been recently used in 
optical cavity quantum electrodynamics in its dissipative 
version ' for studying the decoherence of a single spin 
coupled to a spin bath 10,11 or quench dynamics 12 . Note 
also that, in recent years, the entanglement properties 
of its ground state 13 ' 14 ' 15 - 16 ' 17 ' 18 ' 19 - 20 ' 21 ' 22 as well as the 
finite-size behavior 23 ' 24 ' 25 ' 26 have focused much attention 
on this model. 

An exact solution of this model has been 
derived 27 ' 28 - 29 , but it requires the solution of Bethe-like 
equations, which is more costly in terms of compu- 
tational effort than exact diagonalization. Although 
the low-energy physics of the model has been widely 
studied through different approaches (variational 1 ' 30 ' 31 , 
bosonization 23,32 ' 33 , and coherent states' 5 ' 11 ), its 
high-energy properties have only been very recently 
investigated numerically 35 ' 36,3 ' and several interesting 
features have been revealed. More precisely, for special 
values of the energy, the spectrum has been shown to 
display singularities which are reminiscent of the critical 
point responsible for the well-known quantum phase 
transition at zero temperature. 

In a recent Letter , we proposed a theoretical frame- 
work which allows for an exact computation, in the ther- 
modynamic limit, of the LMG model spectrum for the 
whole range of parameters and leads to a precise descrip- 
tion and understanding of these so-called exceptional 
points. The present paper is an extension of that work, 
in which we will detail some of its main results and ex- 



tend the analysis along several directions (relation to the 
semiclassical treatment, first-order finite-size corrections 
and expectation values of observablcs). 

The paper is organized as follows. In Sec. II, we 
introduce the LMG model and the spin-coherent-state 
formalism' 11 , which is the key ingredient of our ap- 
proach. We then derive the classical energy surface 3 ', 
whose extrema lead to a qualitative phase diagram; these 
extrema are related to the density-of-states singulari- 
ties. In a second step, most importantly, we analyze 
this phase diagram quantitatively. In Sec. Ill, we in- 
troduce the Majorana polynomial and map the time- 
independent Schrodinger equation onto a first-order non- 
linear (Riccatti-like) differential equation. In Sec. IV, 
we give solutions of this equation in the thermodynamic 
limit. This leads to simple expressions of the density of 
states in the whole phase diagram. In Sec. V, we go be- 
yond this limit and compute the leading finite-size correc- 
tions to the density of states. Finally, in Sec. VI, we com- 
pute the expectation values (throughout the spectrum) 
of some spin observables, paying particular attention to 
one region for which spectral subtleties prevent us from 
using the Hcllmann-Feynman theorem. Some technical 
details are given in the Appendix. 



II. COHERENT-STATE REPRESENTATION 
AND CLASSICAL ENERGY SURFACE 



A. Lipkin-Meshkov-Glick model 



The LMG model describes a set of N spin-i parti- 
cles mutually interacting through an (anistropic) XY- 
like Hamiltonian and coupled to an external transverse 
magnetic field h. The Hamiltonian of this system can 
be expressed in terms of the total spin operators S a = 



2 



2^=1 a a/^ where the <r Q 's are the Pauli matrices: 



H 



1 

"N 



( lx S 2 x + 7v S 2 y )-hS z . 



(1) 



In the following, for simplicity, we only consider the 
maximum spin sector s = N/2 with N even. Given the 
symmetry of the spectrum of H, we focus on the param- 
eter range h > 0; [-y^, | ^ j x . Note also that [H, S 2 ] =0 
and [H, e i7T ( s *~ s >] = (spin-flip symmetry). In the stan- 
dard eigenbasis {|s, m)} of S 2 and S z , this latter symme- 
try implies that odd- and even-m states decouple. In the 
thermodynamic limit, both subspaces are isospectral so 
that we further limit the following analysis to the (s + 1)- 
dimensional sector with m even. It is known that H ex- 
hibits a quantum phase transition for h = j x or h = j y . 



B. Coherent-state representation of the spin 
operators 

To determine the spectrum of the Hamiltonian if, it is 
convenient to use a spin-coherent-state representation . 
Let us denote by {|s,m)} the standard eigenbasis of 
{S 2 ,5k} with eigenvalues s(s + 1) and m, respectively. 
The unnormalized spin coherent state |a) is then defined 
as 



e aS +| S ,- S ). 



The scalar product of two such states is 

(l + aa') 2s , 



(a a) 



(2) 



(3) 



where a is the complex conjugate of a. These coherent 
states obey the following closure relation: 



dada (2s + 1) \a)(a\ 
7T (1 + aa) 2 (a\a) 



= 1, 



(4) 



where J dada = JdRe(a) dlm(a). In this representa- 
tion, a quantum state \& (a) = (a| l I / ) is a polynomial in 
a, and the action of the spin operators on ^ translates 
into differential operators: 



S+ = 2sa — a 2 d a , 

S- = d a , 

S z = —s + ad a , 



(5) 
(6) 
(7) 



where S± = S x ± iS y . We shall discuss below the rep- 
resentation of ^(a) in terms of its zeros (the Majorana 
representation) . 



in particular, the location of the quantum phase tran- 
sition. The latter can obtained from an analysis of the 
minima of the variational energy Hq : 



„ , .. Ka\H\a) 
H (a,a) = lim — - 



s^oo s (a a) 



(8) 



2 (l - a 2 a 2 ) h — (a + a) 2 j x + (a — a) 2 7 y 



2(1 + aa) 



(9) 



Note that, in this limit, a classical spin description is 
valid, such that the correspondence between a state |a) 
and a classical vector is simply obtained via a stereo- 
graphic map from the complex plane onto the S 2 sphere 
[with a = e 10 tan(0/2)], leading to the parametrization 

N 

S = — (sin 9 cos 0, sin 9 sin 0, cos 9). (10) 

Here wc shall first be interested in the geometrical 
properties of the whole classical energy surface Ho (a, a). 
Its extrema, obtained by imposing 9 a Ho = <9 Q Ho = 0, are 
given in Table I together with the corresponding energy. 
When one further imposes that a and a be complex con- 
jugate, the configuration space (spanned by the Hamil- 
tonian parameters) is split into distinct regions charac- 
terized by the number of extrema and saddle points in 
Ho (a, a). 

This phase diagram coincides with that derived from 
the analysis of density of states singularities, as done in 
the next section. We shall describe below how far the 
classical analysis can help in understanding the spectral 
results. Note that a related analysis of the classical en- 
ergy surface, including comparisons to numerically de- 
rived spectra, has already been proposed by Castahos 
et al. in terms of the (9, <fi) angles instead of the present 
(a, a). 
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TABLE I: Extrema of the energy surface Hq. 



C. Classical energy surface 

In the thermodynamic limit, a variational description 
of the ground state 1 ' 30,31 , built with respect to the |a) 
states, leads to the dominant behavior of the model and, 



D. Classical description of the phase diagram 

The zero-temperature phase diagram of the LMG 
model is usually discussed in terms of its ground- 
state properties. In this case, only two phases are 
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FIG. 1: Phase diagram in the (7^, jy) plane at fixed h > 
and typical density of states for (7x,7j/,?i) equal to I: (1/2, 
1/3, 1), II: (2, 1/2, 1), III: (5, -3, 1), and IV: (5, 3, 1). 



distinguished 1 ' 24 ' 31 . For h > 7 X (symmetric phase), the 
ground state is unique and lim^oo (S z ) / s — 1, whereas 
for h < 7 X (broken phase), the ground state is twofold 
degenerate and lim^oo (S z ) /s = h/j x . Note that the 
degeneracy in the broken phase arises only in the ther- 
modynamic limit, where the gap between the ground and 
first excited states vanishes exponentially with s. The 
quantum phase transition at h = j x is of second order 
and characterized by mean-field critical exponents as 
well as nontrivial finite-size scaling behavior 23 ' > 25 . 

We have shown " s that, when considering the full spec- 
trum, four different zones arise instead of two, corre- 
sponding to a splitting of the broken phase region into 
three distinct parts characterized by different singulari- 
ties in the density of states (see Fig. 1). 

Note that such singularities have already been pointed 
out in the numerical study of the special case ~f x = 
— 7j/ 35 ' 36 and were called "exceptional points." We em- 
phasize in the present study that these exceptional points 
are associated with saddle points of the energy surface. 
Of course, the absolute minimum (maximum) gives the 
lower (upper) bound of the spectrum. Note that these 
bounds may be degenerate. 

In the thermodynamic limit, to a given energy in the 
spectrum corresponds a level set on Ho (a, a). At that 
energy, the Husimi function local maxima (defined in the 
next section) are known to concentrate along this level 
set, which forms the classical orbit. Singularities of the 
surface (maxima, minima, or saddle points) translate into 
singularities of the level sets (a main ingredient in Morse 
surface theory). This, in turn, affects the density of states 
computation, as illustrated in the next section, and ex- 
plains why the singularities in the Ho (a, a) surface and 
in the density of states are in close correspondance. 

As an illustration, we display in Fig. 2 the classical 
energy surface for (7^ = 5,7 y = 3, h = —1), which is 
precisely the point of zone IV whose density of states is 
shown in Fig. 1. As can be seen, the density of states 
contains two different types of singular points, being the 




FIG. 2: Typical classical energy surface in zone IV (7^ = 
5,7y = 3, h = —1), containing several critical points: two 
minimal points (m)\ two saddle points (S); one local maxi- 
mum (M). It also contains a global maximum, outside the 
range of this plot. The level curves of Ho (classical trajecto- 
ries) are plotted in blue. 



locus of either a divergence or discontinuity. The analysis 
of the classicalenergy surface allows one to qualitatively 
understand all these features. Indeed, it contains two 
absolute minima (denoted m) which provide the lower 
bound of the spectrum (twofold-degenerate ground-state 
energy); two saddle points (denoted s) corresponding to 
the singular behavior of density of states; one local maxi- 
mum (denoted M) which is associated with the disconti- 
nuity, and one absolute maximum, not shown here, giving 
the upper bound of the spectrum. 

The same geometrical analysis can be performed 
throughout the configuration space. A typical classical 
surface in zone I displays one minimum and one max- 
imum, which, respectively, signal the lower and upper 
edges of the spectrum. A zone-II surface has two abso- 
lute minima (corresponding to the broken-phase degen- 
erate ground states), a saddle point (corresponding to 
the density-of-states singularity) , and one maximum (the 
upper spectrum edge) . Finally, a generic zone-Ill surface 
has (again) two absolute minima, two saddle points (cor- 
responding to the two singularities in the spectrum, aris- 
ing at different energies), and two absolute maxima (cor- 
responding to a degenerate upper state). Note that, when 
displayed on the sphere, one recovers the standard result 
for surfaces singularities, which states that the number 
of maxima plus the number of minima minus the number 
of saddle points equals the genus of the sphere, i. e., 2. 

Thus, the analysis of the classical energy surface allows 
us to qualitatively describe the phase diagram shown in 
Fig. 1. However, it does not give any quantitative in- 
formation concerning the density of states. The aim of 
what follows is to develop a reliable method to exactly 
compute the full spectrum of the LMG model. 
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III. MAJORANA REPRESENTATION AND 
SPECTRUM 

A. Majorana polynomial and Majorana sphere 

The first step consists in analysing the eigenstates in 
the spin-coherent-state formalism. Any \^) can be rep- 
resented by its Majorana polynomial ! defined as 



*(a) 




(11) 



(2s)! 



(s — m)\(m + s)! 

CjJ («-«*). 



(s,m\y)a m+s (12) 



(13) 



k=i 



where d ^ 2s is the degree of this polynomial in a (d = 2s 
for a generic state). The roots ctfc of ^>(a) fully charac- 
terize the normalized quantum state |\E') up to a global 
phase. 

It may be more convenient to represent such a state 
|\&) on the so-called Majorana sphere, which can be seen 
as a generalization of the celebrated Bloch sphere used 
for spin-i states. To do so, one first complements the d 
roots vE'(a) with 2s — d roots at infinity in the complex 
plane. Next, the resulting set of 2s complex numbers 
is mapped onto 2s points on the unit sphere by an inverse 
stereographic map. For instance, the basis states |s,m) 
are represented by s — m points on the north pole and 
s + m points on the south pole. Less trivial examples can 
be found in Fig. 3 for eigenstates of H in the zone III. 

Let us also introduce G{a), the logarithmic derivative 

of 



2s f—' a — ah 

fe=i 



(14) 



The 1 /2s factor is here to ensure that G is well behaved at 
the (infinite-s) thermodynamic limit. Let us also define 
the Husimi function associated with a general state \&(a:), 



W\s,(a, a) 



(a|*)(*|a) 



a a 



3 2«[/ Q G(a')da'+ J & G(a' )da' -log( 1+aa)] 



(15) 
(16) 



We shall further need to locate the Husimi function ex- 
trema, which are easily found to satisfy 



G(a) = 



1 + aa 



(17) 



As explained above, for the Hamiltonian eigenstates, 
these maxima converge at the thermodynamic limit, to- 
ward the semiclassical orbits, which are the level sets of 
the classical energy surface Hq. 



B. From Schrodinger to Riccati 

Let us now write the time-independent Schrodinger 
equation H\ty) = E\^S) in the coherent-state representa- 
tion. Using relations (5), (6), and (7), one transforms the 
Schrodinger equation into the linear differential equation 



d a + P (a) 



, Pi(<*) 
(2s) 2 " a 2s 

where e = E/s and 

Po( a ) = ^ a 2 {2s-l)(-fy-j x )-j x -jy 



*(<*)= etf(a), (18) 



Pi(<*) 
P 2 (a) 



a 



2s- 1 
2s 



0?{rix -ly) ~lx ~ly 



(a 2 -l) 2 lx - (a 2 + lf ly 



\-h, (19) 
2/i|,(20) 
(21) 



The next step consists in converting the linear second- 
order differential equation (18) for ^> into a nonlinear 
first-order differential equation for its logarithmic deriva- 
tive G(a), which satisfies the following Riccati- like equa- 
tion 



P 2 (a) 



G'(a) 
2s 



+ G 2 {a) 



+ P 1 {a)G(a) + P (a) = e. (22) 



C. Density of states and poles of G 

The density of states is then obtained from the analysis 
of the poles of the function G. To illustrate the poles 
location, several typical states are displayed in Fig. 3 on 
the Majorana sphere. Each dot represents one pole of G, 
i. e., one Majorana zero a^, which is mapped from the 
complex plane to the sphere by an inverse stereographic 
projection. 

The cornerstone of this study is that, for the LMG 
model, the a^'s spread over two curves Co and C\ in the 
complex plane. In addition, the nth excited state of H 
has 2n poles on C\ and 2(s — n) on Co (thus defining 
both curves). This remarkable property stems mainly 
from existing maps (which may differ between parameter 
space regions) between the LMG model and the problem 
of a particle in an effective one-dimensional potential (see 
the Appendix) . In the latter case, the oscillation theorem 
indexes the excited states by the number of wave- function 
nodes on the real axis. This leads here to (at least one 
set of) zeros lying on simple lines in the complex plane, 
where the pole density varies monotoneously with energy. 

Let us consider the normalized integrated density of 
states, Af(e) <G [0, 1]. We shall enumerate by n the eigen- 
states of increasing energy, starting from n = for the 
ground state to n = s for the highest-energy state. The 
special location of the G poles leads to a simple relation 
between Af(e) and p, the number of poles lying in Ci, 
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FIG. 3: Upper part: representation of the poles of G on 
the Majorana sphere (blue dots) for three typical eigenstates 
computed for h = 1, y x = 5, j y = —3 and s = 20 (zone III 
in Fig. 1). Black lines correspond to the Go branch cuts Co 
and Ci ; orange lines correspond to the classical orbits. Lower 
part: numerical (black dots s — 20) versus analytical (red line 
s — oo) integrated density of states. The two crosses indicate 
the singularities of the density of states NQ U ( — h) and N™(h) 
[Eqs. (39) and (44), respectively] in the thermodynamic limit. 



which reads 
Af(s) 



n 



1 



s + 1 
1 

s+1 



s 

2i7T 



<p_ G(a) da 

JCi 



(23) 
(24) 



where Ci is a contour that surrounds Ci and oriented such 
that M > 0. For the sake of simplicity, we shall further 
consider the density of poles in C\ , called le [0,1], which 
simply reads 



1(e) = 



P_ 

2s 



^— <f> G(a) da. 



(25) 



In general, Eqs. (25) and (22) cannot be exactly solved for 
arbitrary s. The main goal of this paper is to solve these 
in the thermodynamic limit (s — * oo) and to capture the 
leading finite-size corrections in a 1/s expansion. 



IV. THERMODYNAMIC LIMIT 



1. Leading- order expansion for G 



At leading order (l/s)°, Eq. (22) becomes a second-order 
polynomial equation for Go whose solutions are 



Go* (a) 



where 



\ly ~ Jx) +lx+ ly + 2h] ± y/2Q(a) 



2P 2 (a) 



Q(a) = k (a 2 - r£) (a 2 - r\) 

K = ~ {lx ~ ly) + £o) , 



(27) 

(28) 
(29) 



= (-k)- 1/2 ^Jh 2 + lx ly + (7, + ly ) so ± A, (30) 
^{V + 72 + 2 7x e ) (h 2 + 7 2 + 2 ly e Q ) . (31) 



A 



The four roots of Q, ±r±, are branch points of Go- 
The integrated density of states in the thermodynamic 
limit, Ao(eo), now reads 

Afo(e ) = lim N{e) = lim 1(e) = 2 (e a ) (32) 

S— >CX3 s— >oo 

= J_ f da [G+(a)-Ga(a)}. (33) 

zm JCi 

A natural choice for the Go branch cuts is given by 
the curves Co and Ci, on which the G poles accumulate 
as s increases. It indeed corresponds to the direction, in 
the complex plane, for which the quantity computed in 
Eq. (32) is real at each (infinitesimal) step of the inte- 
gration. This latter condition was in fact implemented 
to draw the curves Co and C\ in the different figures. 

In the next section, we analyze in detail the four above- 
mentioned different regions in the phase diagram, in 
terms of Ao(eo), its derivative, and the density of states 
Po(eo) = 9 £o J\fo(eo)- These quantities are, in most cases, 
computed as indicated in Eq. (32). It may happen, as 
noted below, that the C\ curve has a complex shape, 
while Co is simple. Since the integral over all branch 
cuts, corresponding to Co and to Ci, sums to unity, we 
can safely consider the integral over Co, instead of the 
nontrivial one over Ci, and write A/"o(eo) as one minus 
this integral. We also face the case of state degeneracies, 
with corresponding symmetric or nonsymmetric classi- 
cal orbits. Each such orbit is considered separately by 
imposing the analyticity of Go in the region containing 
this orbit, bounded eventually by a closed branch cut on 
the sphere. The related ^(a) is zero along this line and 
can be considered as vanishing outside the considered re- 
gion. This corresponds quite well to the (numerically 
derived) eigenstate in the nonsymmetric case. However, 
in the symmetric case this description fails to reproduce 
the exact eigenstates since the latter is generically a lin- 
ear combination of states located close to the classical 
orbits. 



Let us assume that G and e can be expanded in the 
form 



iGN 



E 



(26) 



2. Analytical expressions of the densities of states 

A precise study of the branch cuts Co and C\ allows one 
to distinguish between five different forms of the density 
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of states (labeled (a), (b), (c), (d), and (e) below) that 
can be expressed in terms of a complete elliptic integral 
of the first kind, 



K(m) 



it/2 



(1 -msin 2 6»)~ 1/2 d6», (34) 



an incomplete elliptic integral of the third kind. 



U{n,(j}\m)= (l-nsin 2 ^)- 1 ^ 
Jq 



m sin 



(35) 

and a complete elliptic integral of the third kind, 
n(n\m) = n(n,7r/2|m). 

Depending on the Hamiltonian parameters, we have 
already distinguished between four different zones, fol- 
lowing the classical surface singularities. We will now 
show how these zones are characterized in terms of the 
density-of-states behavior. Indeed, each time a classi- 
cal surface singularity (maximum, minimum, or saddle 
point) is crossed, the level sets (classical orbits or Husimi 
function local maxima) experience topological changes, 
as well as the integration contours, leading to a new ex- 
pression for the integrated density of states. We now de- 
tail these different expressions, by describing each zone. 

• Zone I: |7 y | < -f x < h. 

Within this range of parameters (which coincides to 
the "symmetric phase" discussed in Sec. II D) the spec- 
trum lies in the interval —h^eo^h and the density of 
states is a smooth decreasing function of the energy as 
can be seen in Fig. 1. The distribution of Majorana poly- 
nomial roots in this zone is similar to that displayed in 
Fig. 3(&). In the complex plane, Co and C\ lie in the imag- 
inary and real axes respectively. The integrated density 
of states is given by 



1 



V2 



2 

a 2 _IlUrt - (36) 

v r_j_ / 



with 



a± = h± ^/j x j y 



(37) 



• Zone II: I7J < h < j x . 

In this region, one must distinguish between two cases 



II (a): 



27* 



£q ^ — h. Co coincides with the 



whole imaginary axis while C\ is made of two discon- 
nected segments in the real axis as depicted in Fig. 3(a). 
Here, the integrated density of states reads 



1 



n 



7rr_. v /27 2; 7y 



(38) 



2 

l- r -± 



i- r 4 



1 — 



n 1 



2 



(1 - M r -) 



— II (b): —h ^ eo ^ h. Cq and Ci are the same as in 
zone I and the analytic expression of the density of states 
is given by Eq. (36). 

These two branches (a) and (6) of the density of states 
diverge at £0 = —h. Indeed, the integrated density of 
states can be simplified into the form 



with 



(39) 



a_ tan 



b + (h) 



a + tan 1 



a+ 



b (h) 



b±{h) = ±(y/h^- Vh^)+^{lx-h) (/i- 7 „),(40) 
b (h) = y/Wy + \j{lx - h)(h--y y ), (41) 

and one can check that p l (—h) = d eo J\f n (eo)\-h di- 
verges. One can further extract the leading behavior of 
the density of states near this point to obtain 



lim Po( £ o) 

eo— >— h 



log \e + h\ 



27r x /(7 x -h)(h- -y y ) 



(42) 



• Zone III: h < — j y < "f x . 

In this region, one must distinguish between three cases 



III (a): 



27* 



£0 ^ —h. Cq and Ci are the same 



as in 11(a), and the integrated density of states is given 
by Eq. (38). 

— Ill (&): — h ^ £0 ^ h. Co and C\ are the same as in 
I, and the density of states N^' '(eo) is given in Eq. (36). 



Ill (c): hs^eo^ 



2 7b 



Cq is made of two discon- 



nected segments on the imaginary axis while C\ coincides 
with the whole real axis as depicted on the Majorana 
sphere in Fig. 3(c). The integrated density of states sim- 
ply reads 



A/o (c) (£ ) = l-CM, 



(a), 



(43) 



where is given in Eq. (38). 

In this zone III, the density of states has two singulari- 
ties at £0 = ±/i. The integrated density of states for these 
energies is given by N^ l (-K) = Af Il (—h) [see Eq. (39)] 
and 



a + tan 



a+ 



bo(-h) 



-a_ tan 



b-(-h) 
(44) 

As done in zone II, one can compute the leading be- 
havior of the density of states near these points and one 
gets 



lim ^( £o ) 



log I 



h\ 



2tt y/- (7* + h) {h + ly ) 



(45) 
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For 7 X = — 7 S , the spectrum is symmetric with respect 
to £ — and the above expression gives the exact loca- 
tion, in the thermodynamic limit, of the so-called excep- 
tional point observed in Refs. ' where a more complex 
diverging behavior was conjectured. 

• Zone IV: h < "f v < j x . 

In this zone the density of states presents three dif- 
ferent regions, of type (d), (e), and (b). The curve C\ is 
more complex here, while Co always lies on a straight line 
in the complex plane. This is why we choose to integrate 
around Co instead of C\ . 



- IV (d): 



Co coincides with 



27x ^ u ^ 2 7m 

the whole imaginary axis while C\ has two disconnected 
branches lying symmetrically on the unit circle with re- 
spect to the imaginary axes. We are here facing a case 
where the classical orbits are related by symmetry [see 
Fig. 4(d)]. One finds, for this region, 



1 



u (— r_) u (r_) 



r+) \/-lxly (h 

£ (r-,y) - 



£ 



(46) 



with 



£ (r-,y) 



n 



n 



n 



n 



u (— r_ 



yu (r_ 
u (— r_ 



yu 
u 



yu 
u 



where 

y 

u (r_) 
- IV(e): - 



yu 



r+ 



sin 1 ^/y 



-r_) 
r_) 



, sin 



, sin 1 y/j/ 



r_ + r + 

h 2 +% 



(47) 

(48) 
(49) 



2jy 



^ So —h. This region shows two 
disconnected classical trajectories not related by symme- 
try (see Fig. 4), corresponding to two qualitatively dif- 
ferent kinds of states which alternate in the spectrum. 
Co comprises two disconnected components lying in the 
imaginary axis, while C\ is still complex and, moreover, 
is different for the two kinds of states. One finds 



A/- (e) M = 1 + 



V2 



IT 1 



- ^O^lxlyK — 



n 



2 \ 


- n \^ir* 


.2 \ ■ 









n 



(50) 



For the critical energy, at the boundary between TV(d) 
and IV(e), the integrated density of states simplifies to 



AC - 



with 



h 2 + l 2 y 
2jy 



c(h) = tan 



^o(-h) = 1- 



1 



[a-c(-h) — a + c(h)~\ . 



3/2 



V /(7,-7,)(7 y 2 -^ 2 ) 



(51) 



(52) 



(53) 



In addition, the density-of-states singular behavior is 
not symmetrical and reads 



lim 



log 



Po 



27y 



= 2 



lim 



(d) 
Po ■ 



, (54) 



(55) 



— IV(6): —h ^ £o ^ h. Co is simply connected and 
lies on the imaginary axes. Like in the previous case, C\ 
is nontrivial (sec Fig. 4). Nevertheless, the expression 
found for Ao in this region coincides with that given by 
Eq. (36). 

We now discuss the particular features found in the 
spectral region IV(e). At e = —h, the density of states 
is discontinuous (see Fig. 1), a fact which can be under- 
stood already from the topological analysis of the clas- 
sical surface Ho. Indeed, the transition from zone (e) 
to zone (b) corresponds to leaving a local maximum of 
Ho (see Fig. 2); therefore, a family of classical orbits no 
longer contributes to the density of states. 

In addition, as opposed to all other regions, the en- 
ergy difference between two consecutive levels, A^ = 
_ £?« ) computed for increasing s, docs not con- 
verge towards the analytical result and, actually, does 
not converge at all. In region IV(e), AW spreads over 
two branches (+) and (— ), depending on the parity of 
the i, which oscillate without converging as s increases, 
as can be seen in Fig. 5. In this case, the gap we com- 
pute, in the thermodynamic limit, is actually the average 
gap, namely A (e ) = §[A(+>(e ) + A(-)(e )]. This is 
clearly to be understood in relation to the existence of 
two kinds of states alternating in the spectrum. Indeed, 
when analyzed separately within each set of states (e + 
or e~ ), the computed energy gaps (between levels j and 
j + 2 in the energy spectrum) converge as s — > oo. In 
addition, both such gaps converge to twice the value of 
Ao(eo) (otherwise the two kind of states would not alter- 
nate as observed numerically). The oscillatory behavior 
noted in Fig. 5 signals an energy drift (with s) of one set 
of energy levels with respect to the other. 
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FIG. 4: Roots of the Majorana polynomial (blue dots) (73, = 
10, 7 a = 5, h = 1, and s = 40), classical orbits (orange curves), 
Co and Ci (black curves), for eigenstates (labelled by n) in 
zone IV (d) (n = 15), zone IV(e) [(e~): n = 25, (e+): n = 
26 ] and zone IV(6) (n = 35). In zone IV(e) two kinds of 
states coexist, of type (e~) and (e + ), associated with the two 
classical orbits nonrelated by symmetry that alternate in the 
spectrum. 



25 
20 
15 
10 




5 = 1040 



FIG. 5: Gap between two consecutive levels as a function of 
the energy in region IV for j x = 15, 7 y = 10 and h = 1. In 
the central region, one sees a real lack of convergence toward 
the red line when increasing s, which is the average gap as 
computed in the thermodynamic limit. 



V. FINITE-SIZE CORRECTIONS 

In the previous section, we have analyzed the thermo- 
dynamic limit of the LMG model spectrum by consider- 
ing the leading terms in the expansion (26) [order (1/s) ]. 
We now express the next-order corrections, which have 



already been shown, at least for the ground state, to dis- 
play nontrivial scaling properties For the sake 
of simplicity, we limit the present analysis to the case 
7^ = 1, and j y = 0. 



A. First-order expansion for G 

Identifying terms of order 1/s in Eq. (22), one obtains 
the following form for the first-order term of G: 



Gf(a) = G 1 (a) + Gf(a), 



with 



G 1 (a) 
Gf(a) 



ha [/i 



1 



1] 



2(l-a 2 )Q(a) 
h (a 2 + 1) + 2 (a 2 



± 



1 £1 



2 (a 2 _ 1) ^/2Q(a) 



(56) 

(57) 
(58) 



Gi is thus an analytic function of a with poles at ±r_ 
and ±r + while Gi has the same branch cuts as Go- 1(e) 
reads, recalling Eq. (25) and developing up to first order, 



1(e) 



1 

2i7r Jc 1 



Go (a) da + 
1 



1 1 



s 2m ./ci 



Gi(a) da, (59) 

Zb(e) + -Zi(e), (60) 
s 

where Jq (e) is given in Eq. (32) and where one can rewrite 
2i(e) = I + -L /" da [G+ (a) - G^ (a)] , (61) 

4 Z17r JCi 

the i coming from the integration over the poles. 

For 72, = l,7j, = 0, one has only zones I and II to 
consider, which focuses the analysis on only two energy 
regions. In zones I and 11(6) one obtains 



I (h + 2 £l )K 



2h n r 



4 7T?"+\ 

whereas in region 11(a) one finds 
1 f 2h 



(62) 



K 1 



r_ (r 

r 2 
- ^2 



1 



- rln [l-rl 



1 - 




+ 



(63) 



Now, for all s, we expect that 1(e) = 2o(eo)> which 
implies, at order 1/s, T\(e — £i/s) = T\(e§) = 0. This 
condition allows one to compute the first-order correction 
to the energy, e\, which is displayed in Fig. 6 (lower left) 
and compares nicely with the numerical values, already 
for small values of s (here s = 50). 
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FIG. 6: Comparison between analytical (red line) and nu- 
merical (s = 50 black dots) results for the (zeroth-order) in- 
tegrated density of states No (upper left) and energy gap Ao 
(upper right) and the first-order finite-size corrections to the 
energy si and to the gap (Ai, lower right). 



of the gap in their vicinity; setting 77 = \h + eq\, one gets 

1 



A(e -> -V 



A(e 



27T^(l-h)h f 1 1 

log 77 



+ 



y/(l- h)h sin" 1 (1- 2h) 
V log 2 T) 



27iV(l - h)h 
log 77 

2^/(1 - Kjhsm^VK 



T) log 7/ 



(66) 



(67) 



Note that the leading term is simply the inverse of po, 
which is given in Eq. (42) and vanishes when 77 goes to 
zero. 



VI. OBSERVABLE EXPECTATION VALUES 



B. Energy gaps 

The gap between two successive levels has already been 
discussed above in the zone-IV case. At the thermody- 
namic limit, it generically reads 



In this section, we discuss the expectation values of 
spin observables for generic eigenstates of the LMG 
model. The simplest way to perform such a calculation 
is to use the Hellmann-Feynman theorem, which relates 
these expectation values to the partial derivative of the 
eigenenergies with respect to Hamiltonian parameters. 
For instance 



A ( £o ) = -tt = wt4' ( 64 ) 
po{e) oNo(e ) 

With the analysis done in the previous section, we can 
now compute finite size corrections to the gap. To first 
order, we obtain 



A = A + -A 1 = A (l + -|^\ (65) 
s V sde J 



The above derived values of £\ allow us to get a closed 
form for Ai, which nicely compares to the numerical val- 
ues, as can be seen in Fig. 6 (lower right) for s = 50. 

The Ai correction is singular at the exceptional points, 
which are, as discussed in Sec. IV, located at £0 = —h. 
Note that Leyvraz and Heiss numerically found a log- 
arithmic singularity at the exceptional points 21 . A re- 
lated feature was already observed for the gap between 
the ground state and the first-excited state 23 ' 24 . In the 
latter scaling hypothesis led to a derivation of the 

first-order correction, showing a A" 1 / 3 behavior. Unfor- 
tunately, the scaling hypothesis cannot be used here at 
the exceptional points. We have determined the behavior 



(®\S Z \$) = -d h E , (*|S ,2 |*) = -2s d lx E. (68) 

As an illustration, we compare in Figs. 7 and 8 three 
cases, computed numerically (at finite s) and via the 
Hellmann-Feynman theorem in the thermodynamic limit , 
i. e. replacing E by ssq. As expected, one can see an al- 
most perfect agreement, except for zone IV(e) discussed 
below. 

Let us still make use of the semiclassical analysis 
discussed in previous sections. The expectation value 
for an observable O reads 41 , at leading order, 

(*idi*) 1 ( T 

{0) = tMt = tJ dt <°Wl°l°W>' (69) 

where T is the period of the classical orbit with en- 
ergy e u and a(t) the solution of the classical dynamics 
equation 42 . 

Let us focus on the (S z ) case. In zone I, it is maximal 
for the ground state. Indeed, in that region, H is mini- 
mum for a = 0, where the classical orbit degenerates to a 
single point at which the ground-state amplitude ^(a)! 2 
is concentrated. As a result, although this true ground 
state differs from the simple fully polarized state 21 , (S z ) 
reaches its maximum value s. 

This also occurs in regions II and III, for energies corre- 
sponding to the exceptional points. Here, classical orbits 
display a characteristic "figure-8" shape, with the val- 
ues of a therefore differing from zero. The saturation 
effect results in that case from the fact that the period 
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III 



FIG. 7: Comparison of expectation values of several observ- 
ables obtained from numerical diagonalizations (black dots) 
and from the Hellmann-Feynman theorem in the thermo- 
dynamic limit (red lines). Plot parameters: s = 60, zone 
I : (7* = l/2,7i/ = !/3,/i = 1), zone II : (j x = 2, -y y = 
1/2, h = 1), and zone III : (-j x — 5, j y — —3, h = 1). 



IV 



-2.0 -1.0 



FIG. 8: Same as Fig. 7, for a typical point in zone IV 
(7z = 5,7 H = 3,h — 1) and s = 60. In the central region 
[zone IV(e)], there is a clear discrepancy between the numeri- 
cal values (black dots) and those derived from the Hellmann- 
Feynman theorem (red lines). 



of the orbit diverges, with a vanishingly small classical 
velocity near a = 0, forcing the expression in Eq. (69) to 
saturate. In all cases except zone IV(e), this latter com- 
putation leads to the same result as that simply obtained 
from the Hellmann-Feynman theorem. 

In zone IV(e), the numerically computed expectation 
values alternate along two distinct curves, differing from 
the Hellmann-Feynman result. This corresponds to the 
already discussed existence, for the same cenergy sq, of 
two kinds of classical trajectories nonrelated by symme- 
try (see Fig. 4). For each numerically derived eigenstate, 
the associated ^(a)! 2 concentrates alternatively near 
one of the two classical orbits. Integrating separately 
along each orbit precisely gives the two branches that 
arc observed numerically (Fig. 8), while the Hellmann- 
Feynman computation leads to an averaged value. 



VII. CONCLUSION 

We have studied in detail the full spectrum of the 
Lipkin-Mcshkov-Glick model by means of a coherent- 



states formalism. In a first step, we simply determined 
the main characteristics of the (zero temperature) phase 
diagram by analyzing extrema and saddle points of the 
classical energy surface. This leads us to distinguish be- 
tween four zones in the phase diagram corresponding to 
various patterns of the density of states whereas the usual 
ground-state criterion leads to only two distinct phases. 

In a second step, we analyzed more deeply the nature 
of the eigenstates in terms of their associated Majorana 
polynomial roots. This enabled us to exactly compute 
the integrated density of states in the thermodynamic 
limit as well as the first finite-size corrections. This re- 
markable result mainly stems from the fact that the roots 
of the Majorana polynomial lies on well-defined curves, 
where their density varies monotoneously with the en- 
ergy. We also clarified the nature of the so-called "ex- 
ceptional" points in the spectrum. 

Finally, we addressed the question of computing 
generic observable expectation values, in particular 
when, owing to subtle spectral reasons, the Hellmann- 
Feynman theorem cannot be used. 

In principle, the same type of analysis could be per- 
formed for any spin Hamiltonian expressed in terms 
of single-spin operators (so-called "collective models"). 
Preliminary investigations of such models with cubic or 
quartic interactions are currently under study. Another 
perspective, also presently under investigation, concerns 
the dynamical properties for evolutions under both fixed 
and variable Hamiltonian parameters. 
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APPENDIX A: MAPPING THE LMG MODEL 
ONTO AN EQUIVALENT ONE-DIMENSIONAL 
MODEL 

The density-of-states calculation given in this paper re- 
lies on the fact that the roots of the Majorana polynomial 
lie on well-defined curves in the complex plane. This re- 
sult stems from the well-known wave-function node oscil- 
lation theorem for one-dimensional systems, which arise 
here via a mapping of the LMG model onto the problem 
of a particle in a onc-dimcnsional potential (sec' 1 for a re- 
view), which we summarize here. A one-to-one relation 
exists between the energy spectrum of the spin system 
and the low-lying quantum states of such a particle. 

We aim to rewrite the equation for the eigenstate ^(a) 
as a Schrodingcr equation for a particle moving in a one- 
dimensional potential. The procedure consists in three 
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steps, given first for the case 7 y < 0. 

1. We change H into an equivalent form such that 
the roots of the Majorana polynomials (nodes of 
the wave-function) which are aligned on the C\ 
curve are sent onto the unit circle. This is achieved 
through the following unitary transformation: H = 



I: 7l = 1/2; 7, = -1/3 II; 7, = 2; J, = -1/3 HI; y x = 5;y, = -2 



e 2 • 



2. The unit circle being parametrized by an angle 0, 
we write $(0) = e" ise *(e ie ) for 9 6 [0,2tt[. 

3. Finally, we define a new function <f>(x), which sat- 
isfies a one-dimensional Schodinger equation and 
such that part of its spectrum is put in one-to-one 
correspondance with the original spin spectrum. 
This is achieved by setting $(0) = e^ x ^(j)[x{9)\ 
where f(x) and x(9) are chosen to suppress the 
first-order derivative in the initial Equation (18) 
for \fi(a) and to set the "mass" term equal to s. 
The resulting Schrodinger-like equation for <f>(x), 
describing a particle in a one-dimensional periodic 
potential, reads 



Y s d 2 M + V(x)<f>(x) 



Ecj>(x) 



(Al) 



Following this procedure, one obtains the effective po- 
tential 



V{x) 



1 



2^ y - 2~j x &n{B\^ x /^ y )- 



A 



h(2s + 1) (7a, - 7„) sn(B\~/ x /j y ) ~ 

[h 2 s + (s + 1)7,7,] cn(B\ lx / ly ) 2 Y (A2) 



with 



B 



K 



(A3) 



Note that V is periodic with period L 



The mapping onto a one-dimensional potential and the 
celebrated node oscillation theorem allows one to sort 
the eigenstates of increasing energy according to their 
number of nodes. Clearly, a <fi(x) node leads to a 'J(a) 
node for the corresponding LMG eigenstate. The first 
(2s + 1) eigenstates of this Hamiltonian H correspond to 
the eigenstates of the LMG Hamiltonian with the same 
energy. Note that, since we focus in this paper on the 
(s + l)-dimensional "even-m" sector, this leads eventually 
to a node number inceasing by steps of 2 for each new 
eigenstate. 

Typical potentials arc shown in Fig. 9, with param- 
eters associated with regions I, II and III of the LMG 
phase diagram. The LMG spectrum corresponds to the 
energies lying between the lower (blue) and the upper 
(red) lines. The qualitative differences between the three 
regions appear clearly here. Indeed, in region I the par- 
ticle moves in a single-well potential whereas it is in a 



v»M o. 




i 0.0 0.5 1.0 1.5 



FIG. 9: Effective one-dimensional potential in the thermody- 
namic limit Vao{x) = lim s _ 00 for j y < and h = 1. Blue 
and red lines are respectively the lower and upper bounds of 
the spin system spectrum Eo = — . 



I: 7 I = l/2;7» = l/3 II: 7x = 3;7» = l/5 



IV: 7l , = 5;7»=3 




-4 -2 2 4 6 




.5-1.0-0.5 0.0 0.5 in 1.5 



FIG. 10: Effective one-dimensional potential in the thermo- 
dynamic limit Voo{x) = linis^oo ZiEZ for 7 H > and h = 1. 
Blue and red lines correspond, respectively, to the upper and 
lower bounds of the energy eo = — in the LMG problem. 



double-well potential in region II. In region III, a higher 
"allowed" energy region appears, with the extended (un- 
bounded) states above the potential barrier. Crossing 
the latter corresponds to the upper density-of-states sin- 
gularity discussed in the text. Note, however, that the 
extended or bounded nature of the eigenstates for this 
equivalent one-dimensional system does not have a direct 
translation into the nature of the corresponding eigen- 
states in the LMG problem. 

Similar transformations can be achieved for positive 7 y 



but in this case, one must consider H 



'He- 



Note the occurence of the minus sign which maps the 
high-energy states of the LMG model onto the low-energy 
states of the particle-problem (and reciprocally). Follow- 
ing steps (2) and (3), one obtains the potential 



V(x) 



27 y cn[C |7„/(7j, - -fx)} 2 ~ 2 lx 
h(2s + 1) {j x - y y ) cn[C |7j,/ (i y - 7,)] - 

(h 2 s + (s + 1)7* 7 S ) sn[C7 | 72/ /( 7y - j x )f 



(A4) 



with 



C = \jlx - ly X. 



(A5) 



-K 



Here, V is periodic with period L = —j=^ 

The effective potentials are displayed on Fig. 10 for zones 
I, II and IV, where some care must now be taken for the 
correspondence with the LMG model. The upper levels 
(close to the upper red line) correspond to the lower levels 
in the LMG case. 
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